Stochastic oscillations in models of epidemics on a network of cities 
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We carry out an analytic investigation of stochastic oscillations in a susceptible- infected-recovered 
model of disease spread on a network of n cities. In the model a fraction fjk of individuals from city 
k commute to city j, where they may infect, or be infected by, others. Starting from a continuous 
time Markov description of the model the deterministic equations, which are valid in the limit when 
the population of each city is infinite, are recovered. The stochastic fluctuations about the fixed 
point of these equations are derived by use of the van Kampen system-size expansion. The fixed 
point structure of the deterministic equations is remarkably simple: a unique non-trivial fixed point 
always exists and has the feature that the fraction of susceptible, infected and recovered individuals 
is the same for each city irrespective of its size. We find that the stochastic fluctuations have an 
analogously simple dynamics: all oscillations have a single frequency, equal to that found in the one 
city case. We interpret this phenomenon in terms of the properties of the spectrum of the matrix 
of the linear approximation of the deterministic equations at the fixed point. 
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I. INTRODUCTION 



Two of the ideas that are currently dominating the dis- 
cussion of modeling epidemic spread are those of stochas- 
ticity and network structure [l|-|6| . Deterministic models 
of the Susceptible-Infected-Recovered (SIR) type have 
a long history 0, [1] and have been thoroughly investi- 
gated [1, along with many extensions of the mod- 
els such as age classes or higher-order nonlinear interac- 
tion terms. Although stochasticity, due to random pro- 
cesses at the level of individuals, and networks, either 
between individuals or towns and cities, were recognised 
early on as significant and important factors, the ten- 
dency was to model them through computer simulations. 
This is not surprising: it is rather straightforward to deal 
with stochastic behavior in simulations, and similarly the 
analytic methods available to investigate complex net- 
works, especially adaptive networks, are limited. There 
has also been a tendency towards developing extremely 
detailed agent-based models to study disease spread \U\ - 
ITij ]. which are the antithesis of the simple analytic ap- 
proach based on the original SIR deterministic model. 

In parallel with these developments, however, there 
have been several efforts to extend analytic studies into 
the realm of stochastic and network dynamics. The 
SIR model can be formulated as an individual-based 
model (IBM) which can form a starting point for both 
an analytical treatment, based on the master equation 
(continuous-time Markov chain) [ll| [l6| , and numerical 
simulations, based on the Gillespie algorithm 117]. The 
analytical studies use the system-size expansion of van 
Kampen to reduce the master equation to the set of de- 
terministic equations previously studied, together with a 
set of stochastic differential equations for the deviations 
from the deterministic result. As long as one is not too 
close to the fade-out boundary, there is no need to go be- 



yond next-to-leading order in the expansion parameter, 
l/yN, where N is the number of individuals in the sys- 
tem. This already gives results which are, in general, in 
almost perfect agreement with the results of simulations 

[HO- 

This approach has been used to study the stochastic 
version of the standard SIR model (l9j | ; the Susceptible- 
Exposed-Infected- Recovered (SEIR) model [2fj| . both 
these models with annual forcing [20, Hl[ , staged- models 
22 , I23L the pair-approximation in networked models 
2J, [25(, amongst others. In this paper we extend the 
treatment to a metapopulation model for disease spread, 
which consists of n cities (labeled j = 1, ... ,71), each of 
which contains Nj individuals. A fraction of the popu- 
lation of city k, fjk, commutes to city j and this defines 
the strength of the link from node k to j in the network 
of cities. We will show that the methods used in the 
case of one city carry over to the case where the system 
comprises of a network of cities, and that a surprisingly 
simple set of results can be derived which allow us to 
make quite general predictions for a class of stochastic 
networked models of epidemics. 

The starting point for our analysis is a specification 
of how commuters move between cities in the network. 
As will become clear, the model we arrive at will not 
depend on the details of how and when these exchanges 
take place. We then write down transition rates for the 
usual SIR process, now taking account the city of origin 
of the infector and infected individuals. From the result- 
ing equation we can immediately find the deterministic 
equations corresponding to the stochastic model when 
Nj — > 00 for all j. Deterministic models of this type 
began to be considered long ago [26( and the existence 
and stability properties of the endemic equilibrium were 
studied for different formulations of the coupling between 
the cities and of the disease dynamics |27H29| . Stochastic 
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effects in these systems have aiso been anaiyzed from the 
point of view of the relation between spatial heterogene- 
ity, disease extinction and the threshold for disease onset 

[ii mil- 

Some rather strong and general results on the unique- 
ness and global stability of the fixed points of the deter- 
ministic model are known (33[- We will use these results 
and then go beyond this leading-order analysis to deter- 
mine the linear stochastic corrections that characterize 
the quasi-stationary state of the finite system. As ex- 
pected, the qualitative predictions of the deterministic 
model are shown to be incorrect; instead large stochas- 
tic cycles are found, although their form is much simpler 
than might naively have been expected. We show that 
this is, in part, a reflection of the special nature of the 
fixed points of the deterministic model. 

The outline of the paper is as follows. In Section UT1 we 
describe the basic model and apply it to the case of two 
cities. The generalization to the n-city case in given in 
Section Mil The results for the form of the sustained os- 
cillations are given in Section HVl and we conclude in Sec- 
tion |Vl Two appendices contain technical details which 
are too cumbersome to include in the main text. 



II. TWO-CITY MODEL 

In this section we will formulate the model when there 
are only two cities; the general n city case described in 
Section Hm does not introduce any new points of principle 
and is easily explained once the two city case has been 
understood. 

The SIR model consists of three classes of individuals: 
susceptibles, infected and recovered. The number of in- 
dividuals in the three classes belonging to city j will be 
denoted by Sj , Ij and Rj respectively. We assume that 
births and deaths are coupled at the individual level, so 
that when an individual dies another (susceptible) indi- 
vidual is born. This means that the number of individu- 
als belonging to any one city, Nj, does not change with 
time, and so the number of recovered individuals is not an 
independent variable: Rj — Nj — Sj — Ij, where j = 1,2 

There are four processes in the SIR model which cause 
transitions to a new state: infection, recovery, death of 
an infected individual and death of a recovered individ- 
ual. The death of a susceptible individual does not cause 
a transition, since it is immediately replaced by another, 
newborn, individual which is by definition susceptible. 
Of the four listed processes, the final three only involve 
one individual and so only involve one city. The transi- 
tion rates are [19| : 

(a) Recovery of an infective individual (and creation of 
a recovered individual) 



(b) Death of an infected individual (and birth of a sus- 
ceptible individual): 

T 3 = T(Si + l,h - l,S2,h\Si,h,S 2 ,h) = fih, 

T 4 = T(S U h, S 2 + l,h- l\S 1 ,I 1 ,S 2 ,I 2 ) = lih- (2) 

(c) Death of a recovered individual (and birth of a sus- 
ceptible individual): 

T 5 =T(S 1 + l,h,S2,I 2 \S 1 ,I 1 ,S 2 ,l2)=KNi-S 1 -I 1 ), 
T 6 = T(S 1 ,I 1 ,S 2 + IJ^SuI^S^h) = (m(N 2 -S 2 - I 2 ). 

(3) 

Here 7 and /i are parameters which respectively charac- 
terize the rate of recovery and of birth/death. 

The infection processes introduce the role of the com- 
muters. We let /21 be the fraction of the population 
from city 1 which commutes to city 2, leaving a frac- 
tion (1 — /21) of the population as residents of city 1. 
Similarly, for commuters from city 2, as illustrated in 
Fig. 1. We note that the number of individuals in city 
j is Mj = (1 - f k j)Nj + f jk Nk, where j f k. We will 
not specify the nature of the commute in more detail and 
assume that the fjk are a property of the corresponding 
pair of cities that defines the overall average fraction of 
time that an individual from one city spends in the other 
city. These coefficients will be taken as a coarse-grained 
measure of the demographic coupling between the cities 
that will be applied to all individuals independently of 
disease status and do not discriminate between different 
types of stays with their typical frequencies and dura- 
tions. 

To see the nature of the infective interactions that oc- 
cur, we first fix our attention on those involving suscepti- 
ble individuals from city 1. There are four types of term: 

(i) Infective residents in city 1 infect susceptible resi- 
dents in city 1. 

(ii) Infective commuters from city 2 infect susceptible 
residents in city 1. 

(iii) Infective residents in city 2 infect susceptible com- 
muters from city 1. 

(iv) Infective commuters from city 1 infect susceptible 
commuters from city 1 in city 2. 

The rates for these to occur according to the usual pre- 
scription for the SIR model [l9( are: 

(i) (3 (l-/ 2 i) Sr (l-/ 2 i)Ji/Mi, 

(ii) P (W21) £1/12/2/^1, 



Ti = T(Si,Ii ~ l,S 2 ,l2\Si,h,S 2 ,h) = jh, 

T2=T(S 1 ,I U S2,l2-l\Sl,h,S2,l2)=jl2. (1) 



(iii) /3/2iSi (I-/12U2/M2, 



(iv) /3/21S1/21VM2, 
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FIG. 1: A fraction fjk of residents of city k commute to city 
j, where j, k — 1,2. 



where (3 is the parameter which sets the overall rate of 
infection. Adding these rates together we obtain the total 
transition rate for infection of Si individuals as 
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A similar analysis can be made for the transitions in- 
volving susceptible individuals from city 2. Putting these 
results together we obtain the transition rates for infec- 
tion as 

(d) Infection of a susceptible individual: 

T 7 = T{Si -l,h + 1, S 2 , I 2 \Si,h,S 2 ,I 2 ) 



Sih 



Sih 



P\c n — + c 12 - N2 
T(Si,h,S 2 -l,h + l\S 1 ,h,S 2 ,h) 



C21" 



S_2h 

Ni 
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S 2 I 2 
No 



(4) 



where 



Cll 



C12 



C21 



f'22 



/21 



(I-/21) 2 

1 - /21 + fuq hi + (1 - fi2)q ' 

(1 - f 21) fm /2i(i - /i2)g 

1 - /21 + fuq .hi + (1 - fi2)q ' 

(l -f 12) hi | /l 2 (l-/2l) 

/21 + (1-/12)9 i-/2i+/i2g' 



(1 - /i2) 2 <7 



fhq 



hi + (1-/12)9 I-/21 + /12?' 



(5) 



and q = N 2 /N\ . We assume that iVi and N 2 are not too 
different, so that q is neither very small nor very large. 

The model is defined by the transitions rates in 
Eqs. ([T])-(|3}. It is interesting that the transitions due 
to infection depend on the fractions fjk only through the 



constants Cjk defined in Eq. ([5]). Other ways of account- 
ing for commuting individuals would typically still give 
rise to the form given in Eq. ((4]), but with the constants 
Cjk defined in a different way. 

Since our counting of the ways that infection takes 
place was exhaustive, we expect that the constants Cjk 
are not independent. It is straightforward to check that 
they obey the following relations: 



Cll + C12 = 1, c 2 i 



C22 



1, C12 = C 2 lQ. 



(G) 



So there are only two independent parameters in addi- 
tional to the usual SIR parameters /3,7 and /x found in 
the single city case, and we choose these to be C12 and 
q = N 2 /Ni. Our results will be given in terms of these 
two parameters. It is easy to see that, for each q, the 
range of C12 is the interval [0, q/(q + 1)] where the max- 
imum is attained for hi + /12 = 1- While exploring 
the general behavior of the system we will consider the 
Cjk independently of the underlying microscopic model 
as positive parameters that take values in the wider ad- 
missible range defined by the constraints (|6]). 

Having specified the model it may be investigated in 
two ways as indicated in Section Wl First, it can sim- 
ulated with Gillespie's algorithm [17[, or some equiva- 
lent method. Second, it can be studied analytically by 
constructing the master equation and performing van 
Kampen's system size expansion on this equation. This 
will be the main focus of this paper. For notational 
convenience we will denote the states of the system by 
(7 = {S\, h, S 2 , 1 2 }, recalling that the number of recov- 
ered individuals from each city may be written in terms 
of these variables. The master equation gives the time 
evolution of P(cr, t), the probability distribution for find- 
ing the system in state a at time t. It takes the form 

=EE [Ta(vW)P{v',t) - T a (a'\a)P(<j,t)} , 

a' a—1 

(7) 

where T a (a\a'), a = 1, . . . , 8 are the transition rates from 
the state a' to the state a given explicitly in Eqs. ([TJ-Q. 

The full master equation ([7]) cannot be solved, but the 
van Kampen system-size expansion when taken to next- 
to-leading order usually gives results which are in excel- 
lent agreement with simulations. We will see that this 
will also be the case in the extensions of the method 
which we are exploring in this paper. The system-size 
expansion starts by making the following ansatz [l5| : 



Sj=N jt 



Nj> 2 ^ 



h = N ih 



AT 1 / 2 



(8) 



where j = 1,2. Here Sj — limjv ->oo Sj/Nj and ij — 
lim^v ->oo hlNj are the fraction of individuals from city 
j which are respectively susceptible and infected in the 
deterministic limit. The quantities Xj and yj are the 
stochastic deviations from these deterministic results, 
suitably scaled so that they also become continuous in 
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the limit of large population sizes. The ansatz ((8} is 
substituted into Eq. ([7]) and powers of y/ff] on the left- 
and right-hand sides matched up. The leading order con- 
tribution gives the deterministic equations of the model 
and the next-to-leading order linear stochastic differen- 
tial equations for Xj and yj. We shall not describe the 
method in great detail, since it is described clearly in van 
Kampen's book [l5j and in many papers, including sev- 
eral on the SIR model [1, Ojl, [22| . Instead we will outline 
the main results of the approximation in the remainder of 
this section, and give some explicit intermediate formulas 
in Appendix [A"l 

The deterministic equations which are found to first 
order in the system-size expansion can also be obtained 
by multiplying Eq. ([7]) by Si,h, S2 and 1% in turn and 
then summing over all states a. This yields 



dsi 
~dt 
ds 2 
~~dt 
dii 
~dt 
di 2 
~dt 



-f3si (cnii + c i2 i 2 ) + m(! - si), 
~/3s 2 {c 2 ih + c 22 i 2 ) + m(! - s 2), 
/3si (cnii + c i2 i 2 ) - (7 + fj)h, 
I3s 2 (c 2 iii + c 22 i 2 ) - (7 + M)*2< 



(9) 



For the case of cities with equal population sizes, these 
have been previously found and analyzed in [281 ] . In the 
context of the present work we are mainly interested in 
the fixed points of these equations. We will not discuss 
these here, instead we will wait until Section UITl where 
the case of n cities will be discussed when we can give a 
more general treatment. 

Of more interest to us in this paper are the variables 
Xj and yj which describe the linear fluctuations around 
trajectories of the deterministic set of equations © . For 
convenience we will introduce the vector of these fluctu- 
ations z = (xi,X2, 2/1,2/2)- Our focus will be on fluctua- 
tions in the stationary state, that is, about the non-trivial 
fixed point of the deterministic equations (which we will 
show in the next section is unique). The fluctuations ob- 
tained through the system-size expansion obey a linear 
Fokker-Planck equation, which is equivalent to a set of 
stochastic differential equations of the form [l6[ 



^T= Y.AjKZK + TJjit), J = l,. 



(10) 



K=l 



where rjj(t) are Gaussian noise terms with zero mean 
which satisfy {r]j(t)riK{t')) = BjKS(t — t'). Since the 
fluctuations are about the fixed point, the 4x4 matri- 
ces A and B are independent of time, and completely 
characterize the fluctuations. They are given explicitly 
in Appendix [Al 

The fluctuations will be analyzed in detail in Sec- 
tion IIV1 when they will also be compared to the results 
of numerical simulations. Before discussing this, we will 
generalize the discussion of this section to an arbitrary 
network of n cities. 




FIG. 2: Individuals commute between n cities, illustrated for 
a particular network when n — 4. 



III. ARBITRARY NETWORK STRUCTURE 

In this section we will generalize the content of Sec- 
tion QI] to n cities and also find the fixed points of the 
deterministic dynamics in this case. 

A. n-city model 

We use the same notation as in Section UH labeling the 
cities by j and k which now run from 1 to n. It will be 
convenient to introduce the quantity 



fj ~ fkj ) 



(11) 



so that the number of individuals in city j may be written 
as 



M } = 



N, 



k^j 



(12) 



There are, once again, four types of term in the process 
of infection (see Figure 2) and we again fix our attention 
on those involving susceptible individuals from city 1: 

(i) Infective residents in city 1 infect susceptible resi- 
dents in city 1. This gives a rate of (3(1 — fi)S±(l — 
/i)/i/Mi. 

(ii) Infective commuters from city j, j = 2, . . . , n, infect 
susceptible residents in city 1. This gives a rate, 
summing over all j, of /3(1 - /i)Si Y^j^i fij l j/ M i- 

(iii) Infective residents in city j, j = 2, . . . , n, infect sus- 
ceptible commuters from city 1. This gives a rate, 
summing over all j, of - fiVjffrSt/Mj. 

(iv) Infective commuters from city k (including 
city 1) infect susceptible commuters from city 
1 in city j. This gives a total rate of 
PEfrifjiSiEk&fikh/Mi. 
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Since the transition rates for recovery and birth/death 
are simple extensions of those for two cities we can now 
write down the transition rates for the n-city model as: 

(a) Recovery of an infective individual (and creation of 
a recovered individual) 

Tj = T(Si,h, Sj,I 3 1 S n , I n \a) = jlj, (13) 

(b) Death of an infected individual (and birth of a sus- 
ceptible individual): 

n+j = T(Si,I±, . . . , Sj + 1, Ij — 1, . . . , S n , InW) = 1^1 j, 

(14) 

(c) Death of a recovered individual (and birth of a sus- 
ceptible individual): 

Tin+j = T(Si, Ii, . . . , Sj + 1, Ij . . . , S n , I n \a) 

= n{Nj - Sj - Ij), (15) 

(d) Infection of a susceptible individual: 

T'in+j = T{Sl, I\, ■ ■ ■ , Sj — l,Ij + 1 . . . , S n , I n \(T) 



T, 



= Pj2 c J k 

k=l 



Sjlk 



(16) 



where a = {Si, 1%, . . . , Sj,Ij . . . , S n , I n } and where j = 
1, . . . , n. The coefficients Cjk in Eq. (1161) may be read 
off from the terms (i)-(iv), but they are sufficiently com- 
plicated to write down in full that we only list them in 
Appendix [B] In that Appendix we also show that rela- 
tions between the Cjk, analogous to those given in Eq. (|6]) 
for the two-city case hold, and are given by 

c jk = 1; c kj = (-M Cj k ] j, k = 1, . . . , n. (17) 

So in the n-city model, there are n(n — l)/2 indepen- 
dent coupling parameters Cjk and (n — 1) parameters for 
city sizes in additional to the usual epidemiological pa- 
rameters. Note that if all city sizes are equal the second 
relation in Eq. (|17[) reduces to Ckj = Cjk- This symmetry 
will be used in the subsequent analysis. 

Following the same path as in Section HH having speci- 
fied the model by giving the transition rates, we move on 
to the dynamics. The process is Markovian and so satis- 
fies the master equation ([7]) except now the sum on a goes 
from 1 to An. As detailed in Appendix [A] invoking the 
van Kampen ansatz ([5]) gives the following deterministic 
equations to leading order: 

ds - 

-JT = -Psj^2c jk i k + fi(l- Sj), 



dt 
dt 



k=l 



(18) 



where j = 1, . . . , n. At next-to- leading order the fluctua- 
tions are found to satisfy the linear stochastic differential 
equation (flQ|) . but with J,K = 1, . . . , 2n. The two ma- 
trices A and B are given explicitly in Appendix |XJ They 
are independent of time, since they are evaluated at the 
fixed points of the dynamics (|18l) . For the rest of this sec- 
tion we will investigate the fixed point structure of these 
equations. 



B. The fixed points of the deterministic equations 



The fixed points of the deterministic equations ([T 
will be denoted by asterisks. Adding the two sets of equa- 
tions we immediately see that 



(7 +/*)*; = a* j = i,. 



.,71. 



(19) 



Using this equation to eliminate the i*, and also using 
Eq. (fr7|). one finds that 



(13 + 7 + M) 



fe=l 



(7 + M), 3 = 1 , 



Two fixed points can be found by inspection. 



(20) 
First, 



suppose one of the i* is zero, for instance i| = 0. Then 
from Eq. (Til?j) sj = 1. From Eq. (TI5)) we see immediately 
that J2k=i c iki*k — 0. Since the coefficients Cjk are non- 
negative (see Appendix [Bj , then i£ = for all k as long 
as cik 7^ 0. Using the i£ which are zero as input into 
Eq. (fT5|) , in the same way as we did originally for i\, we 
see that as long as the cities are connected by non-zero 
Cjk, then they will have no infected individuals present. 
From Eq. (TIT))) it follows that s* k = 1 for these cities. 
This is the trivial solution where no infection is present 
anywhere in this cluster of connected cities. We will as- 
sume that all the cities are connected either directly or 
indirectly, so that i* k = Q,s* k = 1 for all k. 

Of more interest is what we will call "the symmetric 
fixed point". This has s k = s*, a constant, for all k. 
From Eq. (fTU|) one sees that the i k are also independent 
of k, and we denote them by i*. Using Eq. (|17p . s* and 
i* are found to satisfy the equations 



s* [(13 + 7 + M ) - Ps*} = ( 7 + M ) . 
(7 + MK =M(l-s*), 



(21) 



which are the fixed point equations for the ordinary 'one- 
city' SIR model @, H- As is well known these may be 
solved to give for the non-trivial fixed point 

* 7 + M .* M [P - (l + M)] / 00 ^ 
s = — - — , % = ~ ; . (22) 



Pin + m) 



k=l 



Due to a remarkable theorem, we can assert that the 
symmetric solution given by Eq. (|22|) is the only non- 
trivial fixed point of the deterministic equations (|18p (33j j . 
This is proved by finding a Liapunov function for the n- 
city SIR model. In fact the result is more general than 
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we require and was proved for the SEIR model; in Ap- 
pendix [B] we give the explicit form of the Liapunov func- 
tion for the SIR model and a brief outline of the proof 
following the argument in Ref. [33[ for this simpler case. 
The theorem also tells us that the non-trivial fixed point 
(f22|) is globally stable. Therefore we can now go on to 
study stochastic fluctuations about this well character- 
ized attractor. 



IV. SPECTRUM OF THE STOCHASTIC 
FLUCTUATIONS 

Based on previous studies of stochastic fluctuations in 
the SIR model in different contexts, we would expect that 
the fixed point behavior predicted in the deterministic 
limit is replaced by large stochastic oscillations [TH, . 
In effect, the noise due to the randomness of the pro- 
cesses in the IBM, sustains the natural tendency for cy- 
cles to occur, and amplifies them through a resonance 
effect. Since the oscillations are stochastic, straightfor- 
ward averaging will simply wipe out the cyclic structure; 
to understand the nature of the fluctuations we need to 
Fourier transform them and then pick out the dominant 
frequencies. 

So we begin by taking the Fourier transform of the lin- 
ear stochastic differential equation Eq. (fT0|) (generalized 
to the case of n cities) to find 

2n 

^2 {-iuSjK - A JK ) z K (u) = fjj(w), J = 1, . ..,2n, 

K=l 

(23) 

where the / denotes the Fourier transform of the function 
/. Defining the matrix — icoSjK — Ajk to be &jk (w), the 
solution to Eq. ([25)1 is 



2n 



K=l 



(24) 



The power spectrum for fluctuations carrying the index 
J is defined by 




0,00^— v 

0.2 0.4 0.6 0.8 1.0 



FIG. 3: (Color online) Power spectrum for the fluctuations 
of infectives from simulation of a three-city model with equal 
population sizes plotted as a function of the frequency v = 
uj/(2n) 1/y. The spectrum shown corresponds to city 1, the 
spectra for the other cities are very similar. Metapopulation 
model parameters: Ni = N2 = N3 = 10 6 , C12 = 0.06, and 
C13 = C23 = 0.02. Epidemiological parameters: 7 = 365/13 
1/y, /1 = 1/50 1/y, and = 17( 7 + fx). 



peak is present for a large range of parameter values. 
An example is shown in Fig. 3, where typical values for 
measles @, [1(1 H3] were chosen for the epidemiological 
parameters (we shall keep these values fixed throughout 
this section). 

To understand how this comes about, we first note that 
the number of peaks in the power spectrum is given by 
the form of the denominator, | det $(w)| 2 ; the numerator 
essentially just shifts the position of these peaks some- 
what. Therefore we can understand the number and na- 
ture of the peaks by studying the eigenvalues of $>jk, 
which are those of the matrix Ajk shifted by — iuj. 

Each pair of complex conjugate eigenvalues of Ajk, 
A c , A*, will give rise to a factor in | det ^(w)! 2 of the form 



2n 2n 

EE* 



(|A C | 2 — uj 2 ) 2 + [2Re(A c )w] , and each real eigenvalue of 
Ajk, X r , yields a factor of the form (A 2 +w 2 ) 2 . Peaks in 
the power spectrum are associated with complex eigen- 
values A c of Ajk with small real parts, and their position 
is approximately given by uj ss Im(A c ). In the trivial case 
of one city, n = 1, Ajk has a pair of complex conju- 
gate eigenvalues A* with Re(A] t ) = — /fyt/(2(7 + fj,)) and 
\X t I = \/ /i(/3 — 7 — fi) (see Appendix IBl. The conditions 



1 



(lo)Bkl ^ or a P ronounce d peak for ui close to lva(\ 1 



(25) 

Since $ = — A, where I is the 2n x 2n unit matrix, 
and since A and B are independent of to, the structure 
of Pj(oj) is that of a polynomial of degree An — 2 divided 
by another polynomial of degree 4n. The explicit form 
of the denominator is | det <J> (oj) | 2 . 

Oscillations with well-defined frequencies should show 
up as peaks in the power-spectrum. The structure of 
the power spectrum described above — with the ratio of 
polynomials of high order potentially giving rise to many 
maxima — might lead us to suppose that the spectrum 
of fluctuations would have a rather complex structure. 
In fact numerical simulations indicate that only a single 



are fulfilled because //, the death-birth rate, is small. 
This carries over to the general n city case since, as 
shown in Appendix [B] A^ always belong to the set of 
eigenvalues of Ajk ■ For the parameter values of Fig. 3 
the numerical values of the common eigenvalue pair are 
Xf = — 0.17 ± i2.99, so we expect a peak to be located 
close tov = Im(A 1 : )/(27r) a 0.48 1/y. 

For large demographic coupling, the n city system will 
behave as well mixed system comprising all the cities and 
we expect to find in that limit a power spectrum similar 
to the case n = 1, where each city contributes propor- 
tionally to its size to the overall spectral density. In the 
opposite limit, the n cities uncouple and we will find for 
each city the power spectrum of the one city case. In or- 
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FIG. 4: (Color online) An Argand diagram of the eigenvalues 
for the two-city model with q = N2/N1 = 3/2 and C12 £ [0, 1]. 
The large black dots are the common eigenvalue pair Xf. The 
sets of smaller dark gray (blue) and light gray (green) dots are 
the remaining eigenvalues computed on a uniform grid of 
values of C12 in the interval. The eigenvalues with Re < —6 
are not shown in the plot, they are found for C12 > 0.15. The 
asterisks show the eigenvalues for the parameter values used 
in Fig. 5. 
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FIG. 5: (Color online) Power spectra for the fluctuations of 
infectives from simulation of the two-city model [(red) dots] 
and analytic calculation (black solid curve) plotted as a func- 
tion of the frequency v = uj/(2n) 1/y. The population sizes 
were chosen to be Ni = 10 6 and N2 = 1.5 x 10 6 so that their 
ratio is 3/2. The coupling coefficient C12 = 0.1. The location 
of the eigenvalues for this choice of parameters is indicated in 
Fig. 4 by asterisks and large dots. 



ImA 



der to understand why additional peaks do not show up 
in simulations for intermediate coupling strengths, it is 
useful to consider the case n = 2, for which the eigenval- 
ues of Ajk can be determined analytically and depend on 
a single coupling parameter c\i and the ratio of the pop- 
ulation sizes q = N 2 /Ni (see Eq. © and Appendix |BJ . 
An Argand diagram of the two pairs of eigenvalues, A^ 
and A2 j for the two-city model is shown in Fig. 4. It 
can be seen that as the coupling increases, follow the 
circle C centered at zero that goes through A^, moving 
away from the imaginary axis. Real and imaginary parts 
become of the same order for very small values of the 
coupling, and so we expect the power spectrum to be 
always dominated by the peak associated with A^ that 
characterizes the spectrum in the uncoupled case. This 
behavior carries over to the n-city case with symmetric 
coupling, for which a complete analysis of the eigenvalues 
of Ajk can also be given, see Appendix IB"1 In particular, 
it can be shown that apart from the common eigenvalue 
pair Xf Ajk, has a single (n — l)-fold degenerate addi- 
tional eigenvalue pair that behaves as a function of the 
coupling parameter as described above for n = 2. 

For the coupling parameter that corresponds to the 
values of A^ marked with asterisks in Fig. 4 and for a 
certain choice of population sizes, the infective fluctua- 
tions power spectra for the two-city model obtained from 
simulations and from Eq. (|2"5j) are shown in Fig. 5. We 
find a nearly perfect match between the results of nu- 
merical simulations and the analytical calculations. In 
agreement with the above argument the power spectra 
of city 1 and city 2 are very similar to the power spec- 
trum of the one city case, which in turn is very similar 
to the spectrum shown in Fig. 3 for 3 cities with small 
coupling. In all cases the functional form of the spectral 
density is dominated by the peak associated with the 
common eigenvalue pair A^. As for the amplitudes of 
the power spectra Pj(y), their ratio with respect to the 
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FIG. 6: (Color online) An Argand diagram of the eigen- 
values for a three-city model with equal population sizes, 
C12 G [0, 0.98] and the other parameters as in Fig. 3. The 
large black dots are the common eigenvalue pair A*. The 
sets of smaller dark gray (blue) and light gray (green) dots 
are the remaining eigenvalues A^ (left panel) and A^ (right 
panel) computed on a uniform grid of values of C12 in the in- 
terval. The eigenvalues with ReA^" < —6 are not shown in 
the plot, they are found for C12 > 0.12. 



one city case, rj(u), decreases as the coupling increases. 
For two cities and q = 1, the power spectra P3 and P4 
of city 1 and city 2 are equal and the relative peak am- 
plitudes r3 i 4(f max ) decrease with the coupling strength 
C12 down to 0.5. For other values of q, as in Fig. 5, the 
different peak amplitudes in two cities reflect the sym- 
metry P 3 (v;a 2l q) = Pi{v\ci 2 /q,l/q). Depending on q, 
the ratio r^^iy) may become even smaller than 0.5, but 
due to the symmetry that relates P3 and P4, the ampli- 
tude of at least one of these peaks is always comparable 
to that of the uncoupled case. More precisely, it is easy 
to check that 1 < r3(v;ci2,q) + r^v; C12, q) < 2, where 
the second inequality is satisfied strictly for C12 = and 
the lower bound corresponds to the large coupling limit 
C12 = 1 and to v = ^ max . 

The general case of three cities with no symmetry can 
also in principle be treated analytically because finding 
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FIG. 7: An Argand diagram of the eigenvalues for a four- 
city model with the coupling strength x £ [0, 0.52]. The large 
black dots are the common eigenvalue pair \ f. The remaining 
eigenvalues A^ 3 4 computed on a uniform grid of values of x 
in the interval are shown as sets of smaller gray dots. As 
in the previous figures we only show eigenvalues whose real 
part is larger than —6. Metapopulation model parameters: 
N2/N1 : N3/N1 : Ni/Ni = 2 : 3 : 4, C12 = 1/y/Jl = 2ci 3 = 

5ci 4 /2 = 5C23/2 = 3c 2 4 = 4C34- 



the eigenvalues of Ajk reduces to finding the roots of a 
fourth order polynomial. However, the problem now de- 
pends on 3 independent coupling parameters and 2 pa- 
rameters for city sizes and closed form expressions are too 
lengthy to be useful. An approximate, concise descrip- 
tion of the behavior of the eigenvalues of Ajk can be 
given in terms of only two parameters that measure cou- 
pling strength and coupling asymmetry, see Appendix IB] 
In this approximation, we assume that all the Cjk, j =/= k, 
arc of order ^//Z an d treat fi as the small parameter of 
the system. Simple expressions for the real parts and 
the absolute values of the additional eigenvalue pairs , 
A^ of Ajk up to terms of order fi can be derived [see 
Eqs. (|B22I) and (|B24j) ]. These show that, in this ap- 
proximation, both eigenvalue pairs behave as described 
for the symmetric case. As the coupling increases, both 
eigenvalue pairs follow the circle C centered at zero that 
goes through A^, moving away from the imaginary axis. 
The real and imaginary parts become of the same order 
within the scope of the approximation. Equation (IE>22|) 
also shows how the asymmetry lifts the degeneracy of 
the two pairs A^, A^. As the coupling increases, the 
two eigenvalue pairs move along the circle C at different 
speeds. We have checked that Eqs. (|B22|) and (|B24j) give 
a good approximation to the exact results in the regime 
when the eigenvalues are complex. 

The same behavior is illustrated in Fig. 6, where a plot 
of the exact solutions for A^ 3 is shown for parameter val- 
ues that correspond to taking those of Fig. 3 and allowing 
one of the coupling coefficients to span the whole admis- 
sible range. One of the eigenvalues is shown only up to 
C12 = 0.12, where its real part becomes smaller than —6. 

In Fig. 7 we show numerical results for the behavior 
of the eigenvalues of Ajk in the case of 4 cities with 
different population sizes and a certain choice of the cou- 
pling coefficients Cjk, j, k = 1, 2, 3, 4. We make use of the 
following notation for the diagonal and off diagonal cou- 



pling coefficients (see Appendix [Bj) : cjj = 1 — Cjj x y/JI 
and Cjk — Cjk x yfii, respectively. We then calculate 
the set of three non-trivial eigenvalue pairs as the cou- 
pling strength x varies in a suitable interval, keeping the 
Cjk fixed. These results suggest that the behavior of the 
eigenvalues of Ajk is essentially given by the description 
of the symmetric case, and that more general couplings 
break the degeneracy as in the case n = 3, with no effects 
in the contributions to the peaks in the power spectrum. 



V. DISCUSSION AND CONCLUSIONS 

In this paper we have extended the analysis of a 
metapopulation model of epidemics into the stochastic 
domain. Frequently epidemic models involving a spatial 
component, such as the interaction between several cities, 
are studied purely deterministically (28l . I29I ] or through 
computer simulations 0, HH, H3] . We have demonstrated 
how a stochastic metapopulation model can be studied 
analytically by using a relatively straightforward exten- 
sion of the methodology which was used to study a well- 
mixed population in a single city. 

We adopted a simple specification of residents and 
commuters in order to set up the model. However, the 
coefficients which appear in the dynamical equations are 
generic and would appear in the same form if residents 
and commuters were included in a different way. It is 
evident that there are many ways of characterizing the 
interchange of individuals between cities which will result 
in the same model; only the identification of the coeffi- 
cients with the underlying structure will be different. 

The deterministic form of the model predicts that the 
system will reach a stable fixed point where the propor- 
tion of infected, susceptible and recovered individuals is 
the same in every city. The stochastic version of the 
model also predicts a clean simple result: that the large 
sustained oscillations which replace the deterministic pre- 
dictions of constant behavior, have a single frequency 
which is the same for every city. Moreover, for small, 
large and intermediate coupling between the cities, the 
form of the power spectrum of these fluctuations is closely 
approximated by the power spectrum of the single city 
system. 

It is remarkable that such a simple result occurs 
in what is a quite complicated stochastic nonlinear 
metapopulation model. We hope to explore the range 
of validity of this result and its robustness to the addi- 
tion of new features to the model in the future. In any 
case, we believe that the work presented here will give a 
firm foundation to possible future work, including com- 
parisons with the data available on childhood diseases. 
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Appendix A: System-size expansion 

Here we give some of the key steps in the application of 
the system-size expansion to the model explored in this 
paper. The method has been extensively discussed in 
the literature [1, [HI H9I - I25I ] , and so we confine ourselves 
to a brief outline and to displaying the most important 
intermediate results in the derivation. We will assume 
that we are carrying out the calculation for the n-city 
case discussed in Section IIIII the corresponding results 
for Section [IT] can be obtained simply by setting n — 2. 

The first point to mention is that there are apparently 
n expansion parameters: {N±, . . . , N n }. The method is 
valid if they are all large and of the same order. More 
formally we can take, for instance, N\ = N as the ex- 
pansion parameter and express all the other Nj in terms 
of it: Nj — Nqj, where the qj = Nj/N, j = 2, . . . ,n are 
of order one. In practice the method seems to work well 
when the Qj are significantly different from one, but this 
has to be checked a posteriori, for instance by comparing 
the analytic results with those obtained using computer 
simulations. In what follows we will not introduce the qj 
explicitly; we will simply take all the iVj's to be of the 
same order in the expansion. 

The van Kampen ansatz Eq. ((5J) replaces the discrete 
stochastic variables a by the continuous stochastic vari- 
ables z and so we write the transformed probability dis- 
tribution P(a,t) as II(z,t). Since this transformation is 
time-dependent, substituting the ansatz into dP/dt on 
the left-hand side of Eq. flTJ gives [H[ 



dP{o,t) _ dn(z,t) 



dt 



dt 



E 



E 

3=1 

<9II(z, t) dij 
dyj dt 



<9II(z, t) dsj 



dt 



(Al) 



The right-hand side of the master equation ([7]) can 
be put into a form from which it is simple to apply the 
expansion procedure. To do this one introduces step- 
operators [15| defined by 

e Sj /('-'Ij •••)<%)•••) 4j Iii • • - > In) 

=f[S%, . . . , Sj i 1, . . . , S n , 1%, . . . , I n ), 
e ^ 1 /(S'i, • ■ • , S n , Ii,...,Ij,...,I n ) 

=f(S 1 ,...,S n ,I U ...,Ij±l,...,I n ), (A2) 



as 



dP(a, t) 



dt 



E 



(>,, - I ) Tj ( ^ I • 



1 ) ?2n+j + f — - — 1 I Ts n+ j 



P(a,t). (A3) 



Within the system-size expansion these operators have 
a simple structure: 



_ ~ NJ P/2 QP _™ 



Nj p/2 dp 



(A4) 



and so all the terms of the right-hand side of Eq. (IA3|) 
may be straightforwardly expanded. Comparing these 
with the left-hand side in Eq. (|A1|) the leading order 
(~ -JNj) yields the deterministic equations given by 
Eq. (|18l) . The next-to- leading order (which is of order 
one) gives a Fokker-Planck equation: 



m=- E Wj [AjKXK ^ + ^ E B 

J,K=l J 



JK 



J,K=1 



a 2 n 

dzjdzR ' 
. (A5) 

The In x In matrices A and B which appear in this 
equation have the following form. Writing A in blocks of 
four n x n submatrices: 



A = 



the elements of these submatrices are 



r a^ 






A^ \ 



(A6) 



A 



jk 

A jk - 



-fiSjh - /3Sjk Cjtit, 



1/2 



SjCjk, 



A3) 
A jk 



= /35jk};Cjtu, 



A 



(4) 



1/2 



SjC jk . (A7) 



Writing B in a similar way to A in Eq. (|A6j) . the elements 
of the submatrices are 



B f k = V s jk(l- Sj)+ pSjk^Sj 



B 



(2) _ R (3) _ 
jk 



B fk = -^jkij - P § jk E s i c J^' 



By k ] = (7 + n) 5, k ij + 138 jk Sjc je i e . 



(A8) 



for a general function / and where j = 1, . . . , n. Using 
these operators the master equation ([7]) may be written 



From Eqs. (IA7j) and (1A8|) it is clear that the matrices 
Aj k and Bj k depend on the solutions of the deterministic 
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equations given in Eq. (|T8|) . However, since we will be 
interested only in fluctuations about the stationary state, 
these matrices are evaluated at the fixed point. Since the 
unique stable fixed point is the symmetric one, the same 
for all cities, the entries (|A7[) and (|A8|) are given by: 



,*(i) 



»,,, /' - ■'< J <>jl, ■ -I,,, ' I S C jk , 



N; 



1/2 



*(3) 



N 



-y* "jic, ^jk ~ ^ \ N k I " w 1 



1/2 



s*c jk - (p + j)Sjk, 
(A9) 



and 



B'W = 2u (1 - s*) S Jk , B*j? = 2 ( 7 + M ) i** Jfc , 



j re 



(3) 



= -i* + /3s*] <5jj 



(A10) 



where we have used the fixed-point equation (f2~Tj) to sim- 
plify some of the entries in Eq. (|A10[) . 

Finally, the Fokker-Planck equation (|A5[) is equivalent 
to the stochastic differential equation (fTU)) . We will work 
with the latter, since we wish to use Fourier analysis to 
analyze the nature of the fluctuations, and since Eq. (fTU)) 
is linear, it can easily be Fourier transformed, as dis- 
cussed in detail in Section ITVl 



Appendix B: Some results for the n-city case 

In this Appendix we give some of the derivations for 
the n-city case discussed in Section IIIII and Section IIVI 
which are too long and cumbersome to be given in the 
main text. 



for j,k = 1, . . . , n and j ^ k. 

To prove the first relation given in Eq. (|17[) , consider 
the sum cjj + J2k^j c i k - The first term in Eq. (|B1| com- 
bines with the first term in Eq. (IB2|) to give (1 — fj). 
The last term in Eq. (IF31|) combines with the last term in 
Eq. ([B2j to give 



ftj fijNj + J2kjt j,i fikNk 



V 

w [(i-iW + E^JWV, 



y : 

[{i-h)N l + Y, m & 



fe m N„ 



= E 



fkj 


Em^fc fkmN m 




[(1- fk)N k + J2 m ^ k fkmN m 



(B3) 



where in the last line we have performed a relabeling. 
Combining the middle term of Eq. (|B2|) with the result 
in Eq. (IB3|) gives 



E 



(i-/fe)^ + E m ^// 



N„ 



% [(i-fk)Nk + E m ^h 



= fj, (B4) 



using Eq. (jlll) . Adding this to the result (1 — fj) found 
earlier proves the result Cjj + E/c^j c jk = 1- 

We also note from Eq. {B2| that Cj k /N k is symmetric 
under the interchange of j and k Therefore 



Cjk 

Nk~ 



2M. 

N, 



(B5) 



which is the second relation in Eq. (|T7|) . 



1. The coefficients c 



2. Uniqueness and stability of the fixed point 



The coefficients Cj k appearing in Eq. (1161) may be read 
off from the four types of term (i)-(iv) given in Sectionflll] 

{1-fjfNj 



[(1 - fj)Nj +Em#i fir- 



E 



Nrr 



flNj 



for j = 1 , . . . , n and 

0--fj)fjkN k 



Cjk 



'^-fj)Nj+Y. m ^jfjmN n 

fkj (1 - fk) N k 



[(l-/fc)^+E m ^/fcr 
ftj fekNk 



£mN m 



(Bl) 



(B2) 



In Section ITU we asserted that the deterministic equa- 
tions (|18j) have a unique non-trivial fixed point, which 
was globally stable. Here we prove this by giving a Lia- 
punov function for the dynamical system in the invariant 
region R = {(si, . . . , s n , h, . . . ,i n ) : < Sj < 1, < ij < 
l,Sj + ij < l,j — l,...,n} where the system is defined. 
This is a modification of the function given in Ref. (33j 
for the SEIR model. The proof assumes that the ma- 
trix of the coupling coefficients Cjk is irreducible, which 
means that any two cities have a direct or indirect in- 
teraction. Otherwise the proof breaks down because the 
n cities may be split into non-interacting subsets, and 
several equilibria may be found by combining disease ex- 
tinction in some subsets with non-trivial equilibrium in 
other subsets. 

Let /3 jk = j3c jk s*i* kl where (s*, . . . , s* , i\, . . . , i*) is a 
fixed point of Eq. (fT8|) , and denote by M the matrix 
defined by M kj = /3 jk J ^ fc, and Efc=i M kj = 0,j = 
1, ...,n. It can be shown ([33j], Lemma 2.1) that the so- 
lution space of the linear equation Mv — is spanned 
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by a single vector (v\, . . . ,v n ), Vj > 0, j 
L(si, . . . , s n , iii ■ • • i in) be defined as 



L = Z^ V 3 ( S J - s *i s i + *j - e j lo § *J ) • 



L has a global minimum in R at the fixed point. Func- 
tions of this form have been used in the literature as 
Liapunov functions for fixed points of ecological and 
epidemiological models, whose variables take only pos- 
itive values (33[. Differentiating L along the solutions 
of E q. dTSI) . and following the proof of Theorem 1.1 in 
Ref. [23|, we obtain 



n. Let Af k ] and AfJ in Eq. ££7) or in A 



L < 



3,k=l 



Sj l k tj 



S j l k % i 



(B6) 



The properties of the coefficients Vj in the definition of 
L play a crucial role in the derivation of the second term 
in this inequality. Use has been made of the identity 

n n n 

X v i Y fakSjik = Y v i(^ + m)*j"> ( B7 ) 

3=1 k=l 3=1 

which in turn uses the fact that Mv = can be written 

as 



3=1 



E 

3=1 



Pcj k Sji%Vj, ,k=l, 



(B8) 



Following Ref. [33], it can then be shown that the right- 
hand side of Eq (IB6|) is strictly negative except at 
(s*, . . . , s* , i*, . . . , i* ). Therefore, L is a Liapunov func- 
tion for this fixed point in R, and the fixed point is unique 
and globally stable. Note that the result also holds when 
the disease transmissibility /3, the recovery rate 7 and 
the birth-death rate fj, are different in different cities, in 
which case the non-trivial equilibrium is in general not 
symmetric. 



3. Nature of the eigenvalues of the matrix A 

In this subsection we will give some results on the 
eigenvalues of A which are required for the discussion 
in Section HVl 

We first recall that A is closely related to the stabil- 
ity matrix of the deterministic equations (fT5|) . In fact, 
in most applications of the system-size expansion they 
are equal. In our case because we have n expansion pa- 



rameters y/Nj, they are not equal, but closely related. 



A simple calculation of the Jacobian, J, from Eq. ()T8|) . 
shows that 

J = 5'" 1 AS', where S 1 = diag (ViVi, . . . , V^V„) . (B9) 

The effect of the transformation is simply that one can 
obtain J from A by omitting the terms (Nj/Nk) 1 ^ 2 in 



4*(2) 



and A 



*(4) 
3 k 



Eq. (|A9|) . This is useful, since it follows from the simi- 
larity transformation (|B9[) that the eigenvalues of A are 
also the eigenvalues of J. So we may study the simpler 
problem of finding the eigenvalues of the Jacobian at the 
symmetric fixed point (|22[) . 

For orientation, let us explicitly calculate the charac- 
teristic polynomial of the Jacobian for the cases of one 
city and two cities. These are 

7; = 1 : Rx(X) = Q-\d 2 \ 2 +diA + d ), (BIO) 

where 

Q = 7 + /j, <2 2 = 7 + n, di = fifi, 

d = M (7 + M )[ / 3-(7 + / x)], (Bll) 

and 

n = 2 : R 2 (X) = Q- 2 {d 2 X 2 + d 1 \ + d^){g 2 \ 2 + 5 iA + .g )- 
Thus, 

R 2 (X) = i?i(A)g- 1 (. 92 A 2 +. 9l A + . go ), (B12) 

where 

92 =7 + M, 9i =/^ + (ci2 + C2i)(7 + m) 2 , 

.go = 0(7 + A0[/8-(l-cia-Qai)(7 + /*)]■ ( B13 ) 

We see that the factor i?i(A) is common, which sug- 
gests that the pair of eigenvalues found in the one 
city case might always be present in the n city case. 
This is easily proved by considering the vector v = 
(vi, . . . , v n , • ■ ■ , v 2n ) T with components satisfying 
Vi = v and Vi +n — v' for i = 1, . . . ,n. Then the eigen- 
vector equation J*v = Av reduces to that for one city as 
required. 

A similar method can be used to find the characteristic 
polynomial for n > 3 cities with equal population sizes, 
where the couplings are equal, that is, 



Cjk = 



1 — (n — l)c, j = k, 
c, j + k, 



(B14) 



where j, k = 1, n. We now take the components of the 
vector to be v± = —v 2 = v, f„+i = —v n+2 = v' , and 
Vi = Vi +n = for i = 3, . . . , n. The eigenvector equation 
J*v = Av now reduces to 



0H 
7 + M 



+ A I v — (1 — nc)(7 + n)v' 



\-^—-fj]v-[ncOy + iJ l )+X]v' = 0. (B15) 
V7 + M / 



Therefore, both solutions of 

Q- 1 (h 2 X 2 + h ln X + h Qn ) = 0, 



(B16) 
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where 



hon = /ti(7 + /i)[/3-(l-nc)(7 + /i)] 



(B17) 



are eigenvalues of J*. This procedure can be repeated 
for 7i—l independent vectors with only four non-zero 
components and the same symmetry as v. Therefore, 
the characteristic polynomial of J*, R n (X), factorizes as 



R n (X) = R 1 (X) [Q- 1 (h 2 X 2 + h ln X + h 0n )Y 



(B18) 



Finally let us consider 3 cities with arbitrary coupling 
and study the eigenvalues of J* in the limit when the off- 
diagonal coefficients Cjk are small and of the same order. 
It will become clear that the coupling range to explore 
corresponds to Cjk of the order of yfjl and it is convenient 
to introduce the notation 

c jk {x) = {\- t » X ^ 3 = k : (B19) 
[c jk x y/Ji, 3TK 

where j,k = 1,2,3 and i is a positive parameter. 
Eq. (|B19j) represents, for each choice of Cjk, a family 
of systems with all the off-diagonal coefficients Cjk of 
the same order, that reaches the zero coupling limit for 
x = 0. The quantity x measures the distance to zero cou- 
pling along each particular family, scaled by y/Jl. Taking 
into account the properties of the matrix Cjk, given by 
Eq. (|17l) . the characteristic polynomial of J* is a poly- 
nomial of degree six that can be expressed in terms of 
this distance x yfji and of three other independent pa- 



rameters. We choose these to be Cjj, j = 1,2,3. We 
know that this characteristic polynomial factorizes as 
i?i(A)(A 4 + p 3 X 3 + p 2 X 2 + p x X + p ), where p 3 , p 2 , Pi 
and po are some coefficients. The roots of i?i(A) are the 
pair of eigenvalues A^ shared by all the characteristic 
equations of this family of systems. The polynomial of 
degree four can be easily found by direct computation. 
For equal city sizes we obtain for the coefficients 



P3 
P2 
Pi 
PO 

where 

(7 = 
P2 = 



7CT x y/Ji+0(n), 

2(/3 - 7) n - 3/4 T 2 p 2 x 2 M + 0(^ 2 ), 

7W-j)a x^ 2 + 0(fi 2 ), 

(/?-7)V + 0(M 5/2 ), (B20) 



Cll + C 2 2 + C33, 



-9 

Cll 



-33 



2(cnC22 + C11C33 + C22C33). 



(B21) 



Keeping only the leading order terms in each of the coef- 
ficients given by Eq. (|B20[) we find a simple approximate 



expressions for the two additional eigenvalue pairs X 2 , 
A 3 . In particular, we find 

Rc(A±) = -Z(v + k) x + 



Re(Xf) 



(cr-fc) x y/ji+0(jj), (B22) 



where 



4(C?! 



P.2 



C11C22 - C11C33 - C22C33). (B23) 



Assuming without loss of generality that £33 > C22 > en, 
k 2 is positive and so k is real. Note that k = in the sym- 
metric case, and in that case (IB22[) coincide in the same 
order of approximation with the roots of Eq. (|B16|) for 
n = 3. The quantities a x y/Ji and k x y/JL that determine, 
in this approximation, the real parts of the two non- 
trivial eigenvalue pairs can be interpreted as the overall 
coupling strength and the coupling asymmetry for a sys- 
tem of family (|B19|) . We also find for the absolute value 
of the eigenvalues 



|A 2 ± 3 l = V / ^7v / M + 0(M) ; 



(B24) 



which shows that, for all families of the form Eq. (|B19[) . 
the eigenvalues A^ 3 of J* move close to the circle C in 
the complex plane centered at zero that goes through 
A^. For arbitrary city sizes, the same calculation can 
be carried out t o find that Eq. (|B22j) and Eq. (|B24j) still 
hold, with Eq. (|B23|) replaced by 



k 2 = a 2 + 1 + g21+g31 fc 2 , 



(B25) 



921(731 

where qjk — Nj/Nk and 

fc 2 = C 2 1 + (c 2 2g21-C33g 3 i) 2 ~2c 11 (c22g21+C3393l)- (B26) 

The behavior of the two non-trivial eigenvalue pairs 
along a family (IB19[) can be described, in this approx- 
imation, in terms of the two parameters a and k that 
characterize the family and of the scaled distance x. As 
x increases away from zero, both eigenvalue pairs move 
along C with speeds whose ratio is given by (a+k) / (a— k). 
The parameter k that measures the asymmetry of the 
coupling causes the splitting of the two pairs with re- 
spect to the degenerate, symmetric case. The first pair 
to reach the real axis does so for x = 4y / /3 — 7/(7(c+A;)), 
which lies within the scope of the approximation. From 
then on the two real eigenvalues keep changing with x in 
such a way that the square root of their product verifies 
the constraint Eq. (IB24[) until for large x the approxima- 
tion breaks down. 
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